Градиентный спуск: \[\theta_{k+1}=\theta_k-\eta\nabla f\left(\theta_{k}\right)\]
Стохастический градиентный спуск: \[\theta_{k+1}=\theta_k-\eta\hat{\nabla f\left(\theta_{k}\right)}\]
Рассмотрим несколько методик осуществления шага и получения ответа
Обновление learning rate:
Для демонстрации SGD в R будем использовать функцию sgd из пакета sgd
sgd(formula, data, model, model.control = list(), sgd.control = list(...))
Возвращаемый объект будет содержать поля
Функция позволяет оценивать параметры следующих моделей:
Для сравнения расширений SGD построим линейную модель \(y=x_1\beta_1+x_2\beta_2+\varepsilon\) и рассмотрим изменение оценок параметров во времени:
set.seed(42)
N <- 10000
sigma <- 0.1
beta <- c(5, 1)
start <- c(-15, -15)
x.1 <- rnorm(N)
x.2 <- rnorm(N)
y <- beta[1] * x.1 + beta[2] * x.2 + rnorm(N) * sigma
df <- data.frame(x.1=x.1, x.2=x.2, y=y)
methods <- c('sgd', 'nesterov', 'ai-sgd', 'implicit', 'momentum')
models <- lapply(methods, function(m)sgd(y~x.1+x.2+0, data=df, model='lm', sgd.control=list(method=m, start=start, shuffle=F, npasses=3, pass=T)))
dd <- do.call(rbind, lapply(1:length(models), function(mi){m<-models[[mi]]; data.frame(xy=rbind(start, t(m$estimates)), method=rep(methods[mi], length(m$pos)+1))}))
colnames(dd) <- c('beta_1', 'beta_2', 'method')
Пути, полученные различными алгоритмами и линии уровня целевой функции
x.min <- -20
x.max <- 20
xy <- meshgrid(seq(x.min, x.max, by=0.1))
X <- as.vector(xy$X)
Y <- as.vector(xy$Y)
Z <- sapply(1:length(X), function(id)sum((X[id]*x.1+Y[id]*x.2-y)^2/length(x.1)))
h<-xyplot(beta_1~beta_2, data=dd, type='l', groups=dd$method, auto.key=T, scales=list(x=list(limits=c(x.min, x.max)), y=list(limits=c(x.min, x.max))), main='SGD paths & target function levels')
h<-h+contourplot(Z~Y*X, region=F,col.regions=brewer.pal(11,'Spectral'), colorkey=F)
h<-h+xyplot(beta[1]~beta[2], pch=4, cex=2, col='black')
h<-h+xyplot(start[1]~start[2], pch=3, cex=2, col='black')
h
coeffs<-lapply(1:length(methods), function(id)list(method=methods[id], coeffs = models[[id]]$coefficients))
coeffs
## [[1]]
## [[1]]$method
## [1] "sgd"
##
## [[1]]$coeffs
## [,1]
## [1,] 5.000435
## [2,] 1.001512
##
##
## [[2]]
## [[2]]$method
## [1] "nesterov"
##
## [[2]]$coeffs
## [,1]
## [1,] 4.999665
## [2,] 1.001213
##
##
## [[3]]
## [[3]]$method
## [1] "ai-sgd"
##
## [[3]]$coeffs
## [,1]
## [1,] 4.998543
## [2,] 1.000185
##
##
## [[4]]
## [[4]]$method
## [1] "implicit"
##
## [[4]]$coeffs
## [,1]
## [1,] 4.999301
## [2,] 1.000442
##
##
## [[5]]
## [[5]]$method
## [1] "momentum"
##
## [[5]]$coeffs
## [,1]
## [1,] 4.999665
## [2,] 1.001213
Для оптимизации произвольных функций методом стохастического градиентного спуска можно воспользоваться обобщённым методом моментов, определив градиент функции;
Например, для функции Розенброка (\(f\left(x, y\right)=\left(1-x\right)^2+100\left(y-x^2\right)^2\)):
set.seed(42)
N <- 10000
start <- c(-1.5, 3)
beta <- c(1, 1)
sigma <- 0.1
x<-matrix(rnorm(N*2), ncol=2)*sigma
y<-as.matrix(rep(NA, N), ncol=1)
gr<-function(theta, xx) {
x<-theta[1]
y<-theta[2]
return(as.matrix(c(400*x^3-400*x*y+2*x-2+xx[1], 200*(y-x^2)+xx[2])))
}
methods <- c('sgd')
models <- lapply(methods, function(m)
sgd(x, y=y, model='gmm', model.control=list(gr=gr, nparams=2), sgd.control=list(method=m, start=start, shuffle=F, npasses=6, pass=T, lr='adagrad')))
dd <- do.call(rbind, lapply(1:length(models), function(mi){m<-models[[mi]]; data.frame(xy=rbind(start, t(m$estimates)), method=rep(methods[mi], length(m$pos)+1))}))
colnames(dd) <- c('beta_1', 'beta_2', 'method')
x.min <- -2
x.max <- 4
xy <- meshgrid(seq(x.min, x.max, by=0.1))
X <- as.vector(xy$X)
Y <- as.vector(xy$Y)
Z <- sapply(1:length(X), function(id){x=X[id]; y=Y[id]; (1-x)^2+100*(y-x^2)^2})
h<-xyplot(beta_1~beta_2, data=dd, type='l', groups=dd$method, auto.key=T, scales=list(x=list(limits=c(x.min, x.max)), y=list(limits=c(x.min, x.max))), main='SGD paths & target function levels')
h<-h+contourplot(Z~Y*X, region=F,at=(1:10)*10,col.regions=brewer.pal(11,'Spectral'), colorkey=F)
h<-h+xyplot(beta[1]~beta[2], pch=4, cex=2, col='black')
h<-h+xyplot(start[1]~start[2], pch=3, cex=2, col='black')
h
coeffs<-lapply(1:length(methods), function(id)list(method=methods[id], coeffs = models[[id]]$coefficients))
coeffs
## [[1]]
## [[1]]$method
## [1] "sgd"
##
## [[1]]$coeffs
## [,1]
## [1,] 0.9997515
## [2,] 0.9992494
Сгенерируем новый набор данных с большим количеством параметров и рассмотрим зависимость ошибки от итерации при различном выборе начального learning rate
Сравним методику с адаптивным изменением learning rate и степенным убыванием
set.seed(42)
N <- 10000
p <- 40
sigma <- 0.1
tform <- diag(runif(p, min=-1, max=1)*20) #%*% matrix(runif(p*p), nrow=p)
data <- matrix(rnorm(p * N), nrow = N)
data <- data %*% tform
beta <- rnorm(p)
beta.0 <- rnorm(1)
reference <- c(beta.0, beta)
y <- data %*% beta + beta.0 + rnorm(N) * sigma
data <- data.frame(y=y, x=data)
rates <- 10^((-2:4))
sgd.models <- lapply(rates, function(r)sgd(y~., data=data, model='lm', sgd.control=list( lr.control=c(r, NA, NA, NA), pass=T, npasses=50, start=rep(0, length(beta)+1))))
sgd.models.rmsprop <- lapply(rates, function(r)sgd(y~., data=data, model='lm', sgd.control=list(lr='rmsprop', lr.control=c(r, NA, NA), pass=T, npasses=50, start=rep(0, length(beta)+1))))
Построим зависимость конечной ошибки в оцениваемых параметрах от learning rate
errors <- sapply(sgd.models, function(m)norm(m$coefficients - reference))
errors.rmsprop <- sapply(sgd.models.rmsprop, function(m)norm(m$coefficients - reference))
df <- data.frame(rates=rates, errors.sgd=errors, errors.rmsprop=errors.rmsprop)
df
## rates errors.sgd errors.rmsprop
## 1 1e-02 1.09077874 0.07253628
## 2 1e-01 0.11802718 0.06313638
## 3 1e+00 0.06529551 0.06214990
## 4 1e+01 0.06230022 0.06205905
## 5 1e+02 0.06207420 0.06204986
## 6 1e+03 0.06205138 0.06204894
## 7 1e+04 0.06204910 0.06204885
Можно сделать вывод, что использование специальной методики обновления learning rate позволяет практически исключить зависимость результата от правильного выбора его начального значения.
Ниже приведены зависимость ошибки оценки на очередном шаге для степенного убывания learning rate:
lines <- lapply(sgd.models, function(m)data.frame(pos=m$pos, error=apply(m$estimates, 2, FUN=function(r)norm(as.matrix(r - reference)))))
lines.lr <- lapply(1:length(rates), function(i){l <- lines[[i]]; lr <- rep(rates[i], dim(l)[1]); data.frame(l, learning.rate=lr)})
lines.all <- do.call(rbind, lines.lr)
xyplot(error~pos, lines.all, groups=lines.all$learning.rate, type='l', auto.key=T, scales=list(y=list(log=T)), main='Error vs iteration number: decreasing learning rate')
Ниже приведены зависимость ошибки оценки на очередном шаге для обновления learning rate с помощью rmsprop:
lines.rmsprop <- lapply(sgd.models.rmsprop, function(m)data.frame(pos=m$pos, error=apply(m$estimates, 2, FUN=function(r)norm(as.matrix(r - reference)))))
lines.rmsprop.lr <- lapply(1:length(rates), function(i){l <- lines.rmsprop[[i]]; lr <- rep(rates[i], dim(l)[1]); data.frame(l, learning.rate=lr)})
lines.rmsprop.all <- do.call(rbind, lines.rmsprop.lr)
xyplot(error~pos, lines.rmsprop.all, groups=lines.rmsprop.all$learning.rate, type='l', auto.key=T, scales=list(y=list(log=T)), main='Error vs iteration number: adaptive learning rate')
//: //: Слишком быстрое убывание learning rate может //: {r,warning=F,fig.width=10,fig.height=10} [//]: set.seed(42) [//]: sigma <- 0.1 [//]: x <- seq(-1, 1, by=0.01) [//]: y <- c(abs(x)) [//]: data <- data.frame(x=x, y=y+rnorm(length(x))*sigma) [//]: [//]: rates <- 10^((-5:5)) [//]: [//]: sgd.models <- lapply(rates, function(r)sgd(y~poly(x, 2), data=data, model='lm', sgd.control=list( lr.control=c(r, NA, NA), pass=T, npasses=100, lr='rmsprop'))) [//]: [//]: lm(y~poly(x, 2)) [//]: sgd.models [//]:
Сравним временную сложность SGD с OLS
Так как для задачи с \(N\) наблюдениями и \(p\) параметрами вычислительная сложность имеет порядок \(O\left(Np^2\right)\), для демонстрации превосходства SGD выберем количество параметров достаточно большим.
set.seed(42)
N <- 20000
p <- 1000
sigma <- 1
data <- matrix(rnorm(p * N), nrow = N)
beta <- rep(5, p)
beta.0 <- 5
reference <- c(beta.0, beta)
y <- data %*% beta + beta.0 + rnorm(N) * sigma
data <- data.frame(y=y, x=data)
rates <- 1
time.sgd<-system.time(sgd.models.rmsprop <- sgd(y~., data=data, model='lm', npasses=3, reltol=1e-9, lr='rmsprop', start=rep(0, length(beta)+1), size=100))
time.lm<-system.time(lm.model <- lm(y~., data=data))
residuals.sgd <- as.matrix(data[,-1])%*%sgd.models.rmsprop$coefficients[-1]+sgd.models.rmsprop$coefficients[1]-y
m <- sgd.models.rmsprop
errors.rmsprop <- c(norm(residuals.sgd)/norm(y), time.sgd[3])
m <- lm.model
errors.lm <- c(norm(as.matrix(m$residuals))/norm(y), time.lm[3])
errors <- rbind(errors.rmsprop, errors.lm)
colnames(errors) <- c('error', 'time')
rownames(errors) <- c('SGD', 'OLS')
errors
## error time
## SGD 0.02387713 6.029
## OLS 0.00619758 21.381
Заметим также, что для задачи нелинейной регрессии время, требуемое в SGD на вычисление величины шага не меняется, в то время как нелинейная регрессия перестаёт иметь решение в явном виде (и требует применения итеративной процедуры).
Продемонстриреум отбор признаков с помощью LARS на примере набора данных diabetes (используя пакет lars)
Датасет содержит в себе 10 оригинальных признаков и набор взаимодействий между ними.
data <- as.data.frame(diabetes)
corrplot(cor(data))
Совершим регрессию с помощью LARS:
lars(x, y, type = c("lasso", "lar", "forward.stagewise", "stepwise", trace = FALSE, normalize = TRUE, intercept = TRUE, Gram, eps = .Machine$double.eps, max.steps, use.Gram = TRUE)
Для того, чтобы выбрать число признаков, воспользуемся кросс-валидацией
cv.lars(x, y, K = 10, index, trace = FALSE, plot.it = TRUE, se = TRUE, type = c("lasso", "lar", "forward.stagewise", "stepwise"), mode = c("fraction", "step"), ...)
X <- as.matrix(data[,colnames(data)!='y'])
cv.result <- cv.lars(X, data$y, index=1:50, type='lar', mode='step')
Можно видеть, что минимум достигается на количестве элементов
min.step = cv.result$index[which.min(cv.result$cv)]
print(min.step)
## [1] 16
Оценим коэффициенты модели
model.lars <- lars(X, data$y, type='lar', max.steps=50)
print(model.lars$beta[min.step, model.lars$beta[min.step,]!=0])
## x.sex x.bmi x.map x.glu x2.sex x2.hdl
## -94.904160 501.616811 251.792232 18.130314 -17.408382 -187.827872
## x2.ltg x2.age^2 x2.bmi^2 x2.glu^2 x2.age:sex x2.age:map
## 467.778725 7.608734 38.748381 69.809696 107.727556 30.045318
## x2.age:ltg x2.age:glu x2.bmi:map
## 8.613394 11.600284 85.688057
Также можно построить изменение оценок коэффициентов по шагам:
plot(model.lars)